#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Mar 30 22:15:33 2019

@author: wuhuan
"""

import numpy as np  # 矩阵操作
import pandas as pd # SQL数据处理

from sklearn.metrics import r2_score  #评价回归预测模型的性能

import matplotlib.pyplot as plt   #画图
import seaborn as sns

# 图形出现在Notebook里而不是新窗口
#matplotlib inline
# path to where the data lies
#dpath = './data/'
df = pd.read_csv("/Users/wuhuan/Desktop/FE_day.csv")

#通过观察前5行，了解数据每列（特征）的概况
df.head()

# 从原始数据中分离输入特征x和输出y
y = df["cnt"]

X = df.drop(["cnt"], axis = 1)

#特征名称，用于后续显示权重系数对应的特征
feat_names = X.columns

#将数据分割训练数据与测试数据
from sklearn.model_selection import train_test_split

# 随机采样20%的数据构建测试样本，其余作为训练样本
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=33, test_size=0.2)
X_train.shape

# 线性回归
#class sklearn.linear_model.LinearRegression(fit_intercept=True, normalize=False, copy_X=True, n_jobs=1)
from sklearn.linear_model import LinearRegression

# 1.使用默认配置初始化学习器实例
lr = LinearRegression()

# 2.用训练数据训练模型参数
lr.fit(X_train, y_train)

# 3. 用训练好的模型对测试集进行预测
y_test_pred_lr = lr.predict(X_test)
y_train_pred_lr = lr.predict(X_train)


# 看看各特征的权重系数，系数的绝对值大小可视为该特征的重要性
fs = pd.DataFrame({"columns":list(feat_names), "coef":list((lr.coef_.T))})
fs.sort_values(by=['coef'],ascending=False)

## 使用rmse评价模型在测试集和训练集上的性能，并输出评估结果
from sklearn.metrics import mean_squared_error
mean_squared_error(y_test, y_test_pred_lr)#mse
print('The rmse of LinearRegression on test is', np.sqrt(mean_squared_error(y_test, y_test_pred_lr)))
from sklearn.metrics import mean_squared_error
mean_squared_error(y_train, y_train_pred_lr)#mse
print('The rmse of LinearRegression on train is', np.sqrt(mean_squared_error(y_train, y_train_pred_lr)))


#在训练集上观察预测残差的分布，看是否符合模型假设：噪声为0均值的高斯噪声
f, ax = plt.subplots(figsize=(7, 5)) 
f.tight_layout() 
ax.hist(y_train - y_train_pred_lr, bins=40, label='Residuals Linear', color='b', alpha=.5); 
ax.set_title("Histogram of Residuals") 
ax.legend(loc='best');

#还可以观察预测值与真值的散点图
plt.figure(figsize=(8, 6))
plt.scatter(y_train, y_train_pred_lr)
plt.plot([-3, 3], [-3, 3], '--k')   #数据已经标准化，3倍标准差即可
plt.axis('tight')
plt.xlabel('True cnt')
plt.ylabel('Predicted cnt')
plt.tight_layout()

#岭回归／L2正则
#class sklearn.linear_model.RidgeCV(alphas=(0.1, 1.0, 10.0), fit_intercept=True, 
#                                  normalize=False, scoring=None, cv=None, gcv_mode=None, 
#                                  store_cv_values=False)
from sklearn.linear_model import  RidgeCV

#1. 设置超参数（正则参数）范围
alphas = [ 0.01, 0.1, 1, 10,100]
#n_alphas = 20
#alphas = np.logspace(-5,2,n_alphas)

#2. 生成一个RidgeCV实例
ridge = RidgeCV(alphas=alphas, store_cv_values=True)  

#3. 模型训练
ridge.fit(X_train, y_train)    

#4. 预测
y_test_pred_ridge = ridge.predict(X_test)
y_train_pred_ridge = ridge.predict(X_train)

## 使用rmse评价模型在测试集和训练集上的性能，并输出评估结果
from sklearn.metrics import mean_squared_error
mean_squared_error(y_test, y_test_pred_ridge)#mse
print('The rmse of LinearRegression on test is', np.sqrt(mean_squared_error(y_test, y_test_pred_ridge)))
from sklearn.metrics import mean_squared_error
mean_squared_error(y_train, y_train_pred_ridge)#mse
print('The rmse of LinearRegression on train is', np.sqrt(mean_squared_error(y_train, y_train_pred_ridge)))

#还可以观察预测值与真值的散点图
plt.figure(figsize=(8, 6))
plt.scatter(y_train, y_train_pred_ridge)
plt.plot([-3, 3], [-3, 3], '--k')   #数据已经标准化，3倍标准差即可
plt.axis('tight')
plt.xlabel('True cnt')
plt.ylabel('Predicted cnt')
plt.tight_layout()

#### Lasso／L1正则
# class sklearn.linear_model.LassoCV(eps=0.001, n_alphas=100, alphas=None, fit_intercept=True, 
#                                    normalize=False, precompute=’auto’, max_iter=1000, 
#                                    tol=0.0001, copy_X=True, cv=None, verbose=False, n_jobs=1,
#                                    positive=False, random_state=None, selection=’cyclic’)
from sklearn.linear_model import LassoCV

#1. 设置超参数搜索范围
alphas = [ 0.01, 0.1, 1, 10,100]

#2. 生成学习器实例
lasso = LassoCV(alphas=alphas)

#1. 设置超参数搜索范围
#Lasso可以自动确定最大的alpha，所以另一种设置alpha的方式是设置最小的alpha值（eps） 和 超参数的数目（n_alphas），
#然后LassoCV对最小值和最大值之间在log域上均匀取值n_alphas个
# np.logspace(np.log10(alpha_max * eps), np.log10(alpha_max),num=n_alphas)[::-1]

#2 生成LassoCV实例（默认超参数搜索范围）
lasso = LassoCV()  

#3. 训练（内含CV）
lasso.fit(X_train, y_train)  

#4. 测试
y_test_pred_lasso = lasso.predict(X_test)
y_train_pred_lasso = lasso.predict(X_train)

## 使用rmse评价模型在测试集和训练集上的性能，并输出评估结果
from sklearn.metrics import mean_squared_error
mean_squared_error(y_test, y_test_pred_lasso)#mse
print('The rmse of LinearRegression on test is', np.sqrt(mean_squared_error(y_test, y_test_pred_lasso)))
from sklearn.metrics import mean_squared_error
mean_squared_error(y_train, y_train_pred_lasso)#mse
print('The rmse of LinearRegression on train is', np.sqrt(mean_squared_error(y_train, y_train_pred_lasso)))

# 看看各特征的权重系数，系数的绝对值大小可视为该特征的重要性
fs = pd.DataFrame({"columns":list(feat_names), "coef_lr":list((lr.coef_.T)), "coef_ridge":list((ridge.coef_.T)), "coef_lasso":list((lasso.coef_.T))})
fs.sort_values(by=['coef_lr'],ascending=False)

#还可以观察预测值与真值的散点图
plt.figure(figsize=(8, 6))
plt.scatter(y_train, y_train_pred_lasso)
plt.plot([-3, 3], [-3, 3], '--k')   #数据已经标准化，3倍标准差即可
plt.axis('tight')
plt.xlabel('True cnt')
plt.ylabel('Predicted cnt')
plt.tight_layout()


#结果分析：线性回归模型的特征系数偏大，由于加入了正则项，岭回归和lasso模型特征系数表现更好一些。
#在测试集上，线性回归的rmse为815，岭回归为811，lasso回归为801，lasso模型中将某些影响小的特征权重值直接设为0，较之于岭回归，lasso模型对噪声数据的抗干扰性更好。



